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Abstract. We provide a derivation of a more accurate version of the stochastic Gross- 
Pitaevskii equation, as introduced by Gardiner et at Q. The derivation does not rely on 
the concept of local energy and momentum conservation, and is based on a quasi-classical 
Wigner function representation of a "high temperature" master equation for a Bose gas, which 
includes only modes below an energy cutoff Eji that are sufficiently highly occupied (the 
condensate band). The modes above this cutoff (the non-condensate band) are treated as 
being essentially thermalized. The interaction between these two bands, known as growth 
and scattering processes, provide noise and damping terms in the equation of motion for the 
condensate band, which we call the stochastic Gross-Pitaevskii equation. This approach is 
distinguished by the control of the approximations made in its derivation, and by the feasibility 
of its numerical implementation. 



1. Introduction 

Notwithstanding the large body of work done on the kinetics and dynamics of Bose-Einstein 
condensates fj), there is still a need for a method of treating these aspects which is both 
accurately related to fundamental theory and at the same time able to be implemented 
practically and reliably. Problems for which this would be particularly useful are those 
of condensate growth, nucleation of vortices and the treatment of heating by mechanical 
disturbance. 

In this paper we will develop the theoretical basis for a description based on a 
stochastic differential equation for a quasiclassical field, which arises from a Wigner function 
representation of the quantum field. Descriptions of this kind have been presented previously 
by Stoof 13 and ourselves Q). Related phase space methods have been presented, which use 
the Wigner function either explicitly |4||51|6||71 or implicitly |8l|9l[T0|[ni[l2|[ni with random 
initial conditions to simulate quantum noise have also been presented, and all have shown that 
a significant proportion of experimental reality can be reproduced. 

The method used here can be seen as a unification of the ideas of quantum kinetic theory 
as presented in 1141 [Tsl 1161 with those of the finite temperature Gross-Pitaevskii equation, 
as developed by Davis and co-workers Isl l9l 1101 . The main idea is that the higher energy 
modes of a Bose gas are largely thermalized and can be eliminated, to produce a quantum 
mechanical master equation. This is also the central concept in Stoof's work — we compare 
our methodology with his in Sect. l5.4l 

We do the eUmination in two stages: 

a) We first eliminate modes with a wavenumber k such that |k| > A, where 27r/A is of 
order of magnitude of the range of the interatomic potential. Under conditions normally 
met in Bose-Einstein condensates, these modes have no occupation, and the effect is 
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to remove from consideration the high momentum components which occur during an 
actual collision. This leads to a "coarse grained" quantum field theory which contains no 
high momentum components, and which can quite accurately be described by a simple 
Fermi delta function pseudopotential 1171 1191 l20l . For this kind of elimination there 
is no need to use the many-body T-matrix formulation, as there would be if A were 
much smaller, and we were required to eliminate thermally occupied states; neither is it 
necessary to introduce the Huang-Yang pseudopotential II21I involving the derivative of 
a delta function. 

b) We then separate the remaining momentum range into a low momentum component (the 
condensate band), and a high momentum component (the non-condensate band). We 
treat the condensate band fully quantum-mechanically, while the non-condensate band is 
treated as a bath of thermalized atoms, which in this paper and |151 is considered to be 
unchanging or in 1161 is treated by a kind of quantum Boltzmann equation. 

c) The condensate band contains much more than simply the condensate. Typically it spans 
an energy range of the order of magnitude of twice the chemical potential. In this paper, 
our criterion will be that all modes in the condensate band should be have sufficient 
occupation (which we take to be more than 10 atoms) for us to be able to make the 
approximations necessary for a Wigner function stochastic differential equation to be 
valid. 

d) The resulting stochastic differential equation is very similar to Stoof's in general 
appearance, but has extra noise terms, and no reference to the many-body T-matrix or 
self energy functions. These effects are provided by solving the stochastic differential 
equation itself. It also explicitly contains the projector into the condensate band, 
which ensures that the solutions of the stochastic differential equation remain within 
the condensate band. 

The organization of the paper is as follows. We outline our description of the system in Sect.|2] 
and present the derivation of the corresponding master equation in Sect.|3] We then consider 
two simplifications of the master equation in Sect.|4]which we call: 

a) The "quantum optical" master equation, which corresponds to the methodology of our 
earlier Quantum Kinetic Theory papers II4III5l[T6l : 

b) The "high temperature" master equation, whose validity requires fcsT to be significantly 
larger than the particle or quasiparticle energies described. 

In Sect.|5la Fokker-Planck equation is derived from the latter form of the master equation, and 
this is equivalent to a set of stochastic differential equations that correspond to the stochastic 
Gross-Pitaevskii equation of the title. These stochastic differential equations are very similar 
to those proposed by Stoof |3l and ourselves [l], but differ in kind of noise considered, and 
in the explicit implementation of the projection techniques of Davis et al. I8ll9l ll0l . Finally 
we conclude in Sect.|6l with a discussion of the range of applicability of the methodology 
developed, and suggest systems that should be investigated within this framework. 

2. Description of the system 

2.1. The cold-collision Hamiltonian 

A system of Bose atoms interacting via an interatomic potential u(x) is almost universally 
simplified in order to separate the short distance dynamics of pairs of atoms as they proceed 
through scattering events from what one normally considers to be the interesting collective 
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behaviour of the gas itself, which takes place on a relatively slow time scale and over larger 
distance scales. There are a number of ways in which one can proceed to an appropriately 
"coarse grained" description, in which the details of the interaction are replaced by a single 
parameter a, the scattering length for interatomic collisions. All methods, either implicitly 
or explicitly, introduce a momentum scale hA, above which all degrees of freedom are 
eliminated. For the description in terms of the scattering length to be valid, there must be 
essentially no occupation of the states eliminated, and this requires 

keT <$:h^A^/2m. (1) 

If A does not satisfy this criterion, a description in terms of the many-body T-matrix rather 
than simply the scattering length is required, as is done in the approach of Stoof 1*31. 

The criterion that motion of the particles inside the range of the interatomic potential be 
eUminated leads to the requirement 

A < 27r/ro, (2) 

where ro is the effective range of the interatomic potential. These two conditions together 
obviously lead to a condition on the temperature 

keT < hy2mrl (3) 

Under the further condition that 

aA < 1, (4) 

we can describe the dynamics of the modes with momenta below the cutoff HA by a 
field operator tp{x), which contains modes only with momentum below the cutoff, and a 
Hamiltonian 

H = Hsp + Hi, (5) 
in which the single particle Hamiltonian is 

= J d3xV>t(x) (^-|^V2 + l/(x)^ ^(x), (6) 
and the interaction Hamiltonian 

Hi^^ /d3xV'^(x)^t(x)^(x)V'(x). (7) 



2 „ 

Effectively, the elimination of the modes above the cutoff has made the replacement to the 
interatomic potential m(x) ud{x.), with 

u = Anh^a/m, (8) 

and, it must be emphasized, the resulting field theory has the cutoff. A, so that the 
commutation relations of the field operator are 

[^(x),V'^(x')] ^T'aIx-x'). (9) 

The function Va plays the role of a kind of coarse-grained delta function; however, it is also 
a projector into the subspace of non-eliminated modes. Using the commutation relation (|9j 
the Heisenberg equation of motion for the field operator takes the form 

9V^(x) 
dt 

In practical computations a momentum cutoff is selected by the spatial grid used, and on this 
scale Va appears like a delta function. 

When all of the conditions are satisfied, this Hamiltonian is well-defined, and the 
Born series is both convergent and accurate at first order. 



ih 



= y"dV7'A(x-x')|-^VV(x') + '^"(x')+W^^(x')^(x')VXx')|. (10) 
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Figure 1. Schematic view of the condensate band, the non-condensate band, and eliminated 
states for a harmonic trap. 



2.2. Condensate and non-condensate bands 

We now divide the states of the system into a condensate band Rc and non-condensate band 
Rnc^ and perform a corresponding resolution of the field operator in the form 

?/;(x) = 0(x)+^wc(x). (11) 

We will describe Rc fully quantum-mechanically, while Rnc will be taken as being 
essentially thermalized. The cut between Rc and Rnc is set in terms of the single particle 
energy Ej^, which is such that particles with higher energy than this can be considered to be 
fully thermalized with very little effect on their energies from the condensate band. 

We want to make clear at this stage that Er <^ h^A^/2m, and thus that the division 
into condensate and non-condensate bands is quite independent of the cut at the wavenumber 
A treated in the previous section. This is essentially the procedure followed in our work on 
the finite temperature Gross-Pitaevskii equation I8ll9l ll0l with the use of the contact potential 
approximation Uo6{x) in the spatial representation of the equations of motion. However, in 
jSj the FTGPE in basis notation is written in terms of the the full interatomic potential, and 
Sect. 5.2 of |8| illustrates how certain terms can upgraded to a T-matrix description that is well 
approximated by a contact potential. The point of view adopted in this paper is consistent with 
our use of the delta function potential in quantum kinetic theory I14lll5l [l6i. 

For the purposes of this paper we have to take into account two criteria: 

a) The highest energy states of the condensate band should have an occupation of the order 
of magnitude 5-10, so that we can apply Wigner function methods to all modes of the 
condensate band. In practice this means that 

^>5, (12) 

since this corresponds to the high occupation limit of the Bose-Einstein distribution. 

b) The single particle energy levels of the noncondensate band should be essentially those 
of the trapping potential V{x), which is usually a harmonic potential. In II15I we showed 
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that this is true to about 3%, provided that 

ER/^ic > 3. (13) 
The two criteria together lead to the requirement (assuming that the maximum value of 
fj-c ~ m) that 

ksT/fi > 10, (14) 
and this condition is usually satisfied in practice. 

2.2.7. Definition of condensate band field operators To effect the division between the 
two bands in the many-particle Hamiltonian we first note that the effective potential in the 
noncondensate band is considered to be little different from the trap potential. Thus we 
can expand tp{yi) in trap eigenfunctions, and define ipNci^) ^'^ '^hat component which 
is expressible entirely in terms of eigenfunctions with energy eigenvalue larger than Er. This 
means that 0(x) is automatically defined through il It . and that it can be expressed in terms 
of trap eigenfunctions with energy eigenvalue less than Er. 

The division can be expressed in terms of projectors, whose form will be important in 
numerical implementations, as 

Pjvc(x,x')= ^n(x)K„(x'), (15) 

En > £' i?. 

7'c(x,x') = l-7'jvc(x,x'). (16) 
Here the y„(x) are a complete set of trap eigenfunctions. Thus 

^7vc(x) = J <f^VNc{^,^m^) - Pivc{^'(x)}, (17) 
0(x) = y"dV7'c(x,x')^(x') =7'c{^(x)}. (18) 

2.2.2. Separation of condensate and non-condensate parts of the full Hamiltonian The 

Hamiltonian is written in a form which separates it into three components; namely, those 
which act within Rmc only, those which act within Rc only, and those which cause transfers 
of energy or population between Rc and Rnc- Thus we write 

H = Hnc + Ho + Hi^c, (19) 

in which H^^c is the part of H which depends only on i/jnc, Hq is the part which depends 
only on 0, and that part of the interaction Hamiltonian which involves both condensate band 
and non-condensate band operators is called Hj c- Substituting (?!)(x) + ijjNci^) 

the Hamiltonian, we get 

Hi,c = Hfl + Hf}. + Hf'l, (20) 

where the individual terms in -ff/.c are the terms involving operators from both bands, which 
cause transfer of energy and/or particles between Rc and Rnc- We call the parts involving 
one (f> operator 



Hl]c = J ^'"'^ (i^Nc(^)^Nc(^)^Nc{x)^ 0(x) 



(21) 
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The parts involving two if) operators are called 



h\% = 2u / d3xV-]yp(x)V'jvc(x)(0t(x)<^(x)) 



+ 1 / d3xV;]vcWV'jvcW('^W<^W) 



+ ^ / d3x^,^c7(x)^jvc(x)(0t(x)<^t(x)). 




(22) 



The parts involving three ((> operators are called 



h\% = u j d3x(0t(x)0t(x)0(x))^jvc(x) 

+ u f d3xV'jvcW(</'^W'^W<^W)- 



(23) 



2.2.3. Connection to the notation of Davis et al. Here we briefly make a link between 
the notation of this paper, which corresponds largely to that of the Quantum Kinetic 
Theory papers 1141 [Tsl il6J . and that of the finite temperature Gross-Pitaevskii equation 



i) The condensate band Rnc of Quantum Kinetic Theory roughly corresponds to the 
coherent region C of the finite temperature Gross-Pitaevskii equation papers, although 
there is the additional requirement in the latter that the mode occupations must be large. 
However, this is also a necessary condition for the condensate band in this paper. 

ii) The non-condensate band Rnc was called the incoherent region / in finite temperature 
Gross-Pitaevskii equation papers. 

iii) For the field operators we have the correspondences 



3. Derivation of the master equation 

The description assumes that Rnc is at least locally thermalized, and thus the atoms in these 
levels are treated simply as a heat bath and source of atoms for the levels in Rc- This is 
defined by requiring the field operator correlation functions to have a thermal form locally, 
and to have factorization properties like those which pertain in equilibrium. The precise nature 
of these local equilibrium requirements is specified in Sect. 13.2.31 

3.1. Formal derivation of the master equation 

The derivation of the master equation follows a rather standard methodology, formulated in 
1221 . and as previously introduced in 1141 1151 [T6l . We project out the dependence on the 
non-condensate band by defining the condensate density operator as 



papers |il||l[Tol. 



0(x) ^ ■0JVc(x) ^ 7?(x). 

iv) The notations for the projectors are connected by 



(24) 



T'cCx, x') ^ V, VNci^, x') ^ Q. 



(25) 



PC = TrArc(p), 
and a projector V on the space of density operators by 

Vp = PNC ® T^T^Ncip) = PNC PC 



(26) 



(27) 
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and we also use the notation 

Q=l-V. (28) 
The equation of motion for the full density operator p is the von Neumann equation 

P= -^[Hnc + Ho + Hj,p], (29) 
= [Cnc + Co + Ci)p. (30) 
We use the Laplace transform notation for any function f{t) 

/ c-''f{t)dt. (31) 
Jo 

We then use standard methods to write the master equation for the Laplace transform of 

v{t) = Vp{t) as 

si{s) - v{Q) = Cgv{s) + VCiv{s) + VCi[s - Cnc - Cq - QCi]-^QCivis). (32) 

In this form the master equation is basically exact. We shall make the approximation that the 
kernel of the second part, the [ term, can be approximated by keeping only the terms 
which describe the basic Hamiltonians within Rc or Rnc, namely the terms Cnc and Cq. 
We then invert the Laplace transform, and make a Markov approximation to get 

vit) = {Co + VCi)v{t) + \^Ci ^ dT cxp {(£0 + Ci)t] QCi^ v{t). (33) 
There are now a number of different terms to consider 

3.2. Terms in the master equation 

3.2.1. Hamiltonian and forward scattering terms These arise from the term (£0 + 
VCi)v{t), and lead to a Hamiltonian term of the form 

He ^ Hq + i/forward, (34) 

where the forward scattering term is defined by 

-ff forward = 2w J d^X flNC {^)<I)H^)4>{^) ^ (35) 

where the non-condensate band particle density is 

fi7vc(x) = Ttnc |V'jvcW^wc'(x)/Oivc| ■ (36) 

This represents the condensate band Hamiltonian corrected for the effect of the average non- 
condensate density on the condensate. This correction is usually small, but can be included 
expUcitly in our formulation of the master equation. The corresponding master equation term 
is 

PC I Ham = ~ ^ [HC,PC] ■ (37) 
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3.2.2. Interaction between Rq and Rnc We now examine the terms in Hj^\ as defined in 
nil , which contain one or 0^ (x); explicitly, the term in Hamiltonian can be written 

i^j^^ = J d^xZ3{K.)(t>^x)+h.c. (38) 
In this equation we have defined a notation 

Z-six.) = ut(jl^^{x)i'Nc{^H'Nc{^)- (39) 
Substituting into the master equation (I33> . terms arise of which a typical one is of the form 



^ f d^x f d^x' f rfrTr 



{Zs{x)zlix',T)pNc} <P\^m^',T)pc{t). (40) 



Here we have introduced the notation for an arbitrary non-condensate band operator 

Anc{^, t) = c'^"C*/''AArc(x)c-'^'«^*/'*, (41) 

and for an arbitrary condensate band operator 

Bc(x, t) = e'^^*/''Bc(x)e-'^^*/'"'. (42) 
This notation will be used frequently in the remainder of this section. 

3.2.3. Correlation functions of Rnc The terms involving Z3, Z3 are averaged over pnc, 
which is assumed thermalized, and is therefore quantum Gaussian. This means that: 

i) We may make the replacement 

(Z3(x)Z](x', r)) -> 2(7Aivc(x)Vjvc(x'' ^)) (V'ivc(x)V]vc(x'' ^)) (V']vc(x)V'ivc(x', r)). 

(43) 

(Terms involving (?/)]vc'(x)?/'jvc(x')) etc. do arise in principle, but give no contribution to the 
final result because of energy conservation considerations.) 

ii) The time-dependence is needed only for small r, and in this case we can make appropriate 
replacements in terms of the one-particle Wigner function F(u, K) 

<c (u + I) ^.c (u - f .) ) « ^ /^^^ d^K F(u, K)e— (--)^ 

(44) 



(45) 



where 



u = , (46) 

V =x-x', (47) 

;iw(K,u) = ^— + l/(u). (48) 
2m 

Since the range of all the K integrals is restricted to Rnc^ it is implicit that in all integrals 

;iw(K,u) > Er. 

iii) This approximation is valid in the situation where i^(K, u) is a smooth function of its 
arguments, and can be regarded as a local equilibrium assumption for particles moving in a 
potential which is comparatively slowly varying in space. 
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3.3. Growth terms 

Using J43I44I45> we can now write 

(Z3(x)Z3^(x',r)) = ^ JJJd'K,d'K2d'K3[l + F{u,K,)][l + F{u,K2)]Fiu,K3) 

^ g-i(wi+ii;2-'^3)-r-i(Ki+K2-K3)-v (49) 

where 

nuji = nuj{Ki,u). (50) 

3.3.1. Master equation terms for growth We can write 

(/i(x, r) = exp(-iLcT)0(x), (51) 

with 

LcH^)= -[Hc,H^)]/h. (52) 
This means that (^) can be written as 

J d^u J d\ + (u, V, Lc)0 (u - I) } P, (53) 

with 

G(-) (u, V, c.) = - ^^^^ jjj d^Kid^Kad^Kall + i^(u, Ki)][l + F(u, K2)]F(u, K3) 



X 

,2 



(27r)8;i2 



/ |-,l(l^l+li'2-tJ3-l^)T-l(Ki+K2-K3)-V {54) 

J —00 

d^Kid^Kad^Ksil + F(u, Ki)][l + F(u, K2)]F(u, K3) 



X 5{uJi +UJ2-UJZ- tj)c-HKi+K2-K3)-v^ (55) 

(The approximation made in (I55> is to say dr cxp(— ifir) = it5{Q) + iP/f2 w 7r(5(il), 
i.e., to neglect the principal value integral). 
We will also need 

G(+) (u, V, c^) « - -^^^ jjj d3Kid3K2d3K3F(u, Ki)^^(u, K2)[l + ^^(u, K3)] 

X 5{uJi +0J2-^-i- tj)e-i(Ki+K2-K3)-v_ (5g) 

If the noncondensate band is taken to be in thermal equilibrium 

^^(u, K) = ^ , (57) 

cxp {^^^^^) 1 

there is the relation 

G(+)(u,v,w) = e(^-'"'")/'^'«^G(-)(u,v,c^). (58) 

This will still be true even if the equilibrium is merely local, with /i and T depending on the 
position u. 
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Collecting all relevant terms together, we get a master equation term in the form 
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pel 

growth 



V 

"+2 



J £u J d^v [{g(-)(u,v,Lc)</'(u-|)}pc, 
- Jd'nJ d\ [pc {g(-) (u, V, Lc)^^ (u - |) } , (u + ^ 
+ Jd'nJ rfV[{GW(u,v,Lc)0t (^^ - pc , (u + ^ 
:,{GW(u,v,Lc)0(u-^)},0t 



d^u / d- 



(59) 



3.4. Scattering terms 



These come from terms involving two (f) operators; of these terms, only those involving one <j> 
and one 0^ can yield a resonant term, corresponding to scattering of a non-condensate particle 
by the condensate. The effect of the non-resonant terms is neglected. 



(2) 

Thus we arrive at an approximation to Hj in the form of a term 
hP ^2 J d^xZ2ix)U{x) +h.c., 



(60) 



where 



Z2 (x) UlpNC (x) i^lfc (^) ' (61) 
C/(x) =0t(x)0(x). (62) 

The term analogous to J40> in the master equation becomes, using the notation of (I41I42> . 

-^J dTUix)Uix',T) J d^xj d^x'{Z2ix)zl{x',T))pcit). (63) 



3.4.1. Master equation terms for scattering These arise from the term 



dr ( [Z2 (x) C/(x) , [Z2 (x' , r) C/(x' , r) , p] ] ) 



NC ■ 



(64) 



We define 

M(u,v,c.) = ^^^ J d'K, J d3K2F(Ki,u)[l + F(K2,u)K(K^-K^)- 

X S{lui — CJ2 — w), (65) 
and if F^K, u) corresponds to thermal equilibrium, as in (I57> . this satisfies the relation 

M(u, V, uj) = e-^'^^'"''^M{u, V, -ui). (66) 
The terms in the master equation arising from this part become 

d^u f d^v 



pel 



U 



d-^u / d-^v 



u+^),{Af(u,v,Lc)f/(u-^)}p 
•c {a/(u, V, -Lc)U (u - ^) } , [/ (u + ^) 



(67) 
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3.5. Terms involving three <j) operators 

The part of the interaction Hamihonian can be written 

Hf^ = u j d3xV;^c(x)^Hx)0^(x)<?!)(x) +h.c. (68) 
The master equation term this time takes the form 

dr j d^x j dV0^(x',T)0(x',r)0(x',r)(?it(x)0t(x)0(x) 

X (V'Arc(x)V]vc(^'))Pc(i). (69) 
These terms are probably very small, and will not be included in our analysis. 

3.6. Connection with the finite temperature Gross-Pitaevskii equation of Davis etal. 

The master equation terms described above are related to the terms in the finite temperature 
Gross-Pitaevskii equation, (29a-d) of |8j. We have 

Growth terms < — > UoV {(r')t(x)7)(x)?7(x))} , 

Scattering terms < — > UqV |2'(/;(x)(77^(x)7}(x))} , 

Terms with three operators < — > UqP {2|?/)(x)p(7}(x)) + ■i/;(x)^(7}t(x))} 

We point out that the finite temperature Gross-Pitaevskii equation scattering term also includes 
the forward scattering part of the Hamiltonian in \35\ of this paper The anomalous term of 
the FTGPE, analyzed in Sect. 5.2 of |i8J is mainly responsible for the replacement of the 
true interatomic potential by the contact potential in the equations of motion, and hence 
the introduction of the momentum cutoff A. In this paper, as noted in Sect. l2.2l we have 
already performed this spatial coarse-graining in the Hamiltonian with the cutoff A, and so 
any anomalous terms will be very small. 

3.7. The full master equation and its stationary solution 
The full master equation can now be written 

PC = Pclnam + Pclgi-owth + Pclscatt ' (^0) 

where the individual terms are given by (I37> . \59\ and (I67> . The stationary solution of each 
of the terms in the master equation MQ\ . and therefore of the master equation itself, is given 
by the grand canonical form 

This follows: 

i) From the equality of terms in ( I59> like 

{g(-) (u, V, Lc)0 (u - I) } = Ps (u, V, Lc)0 - ^) } > (72) 

which can be derived using using (I58> . (17 1> and the commutation relation 

[iVc,0]=-0, (73) 
between the field operator and the condensate band number operator 

Nc^ I d^x(l)^x)(f>{x). (74) 
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ii) From the equality of terms in (|^} like 

[m (u, V, Lc)C/ (u - ^) } = {m(u, V, -Lc)(7 (u - I) } . (75) 
which can be derived using i66i . i7l\ and the fact that [U{x.), Nc] = 0. 

3.8. The projection into the condensate band 



In exactly the same way as noted in Sect. l2.1l the projection into the condensate band means 
that the operator Lc has a projected expression when it acts on the condensate band operator 

0(x) 

iLc0(x) = Vc [-^^^^i"^) + ^(^) + w'/'^(x)0(x)0(x)| , (76) 
which arises directly from the fact that the field operator commutation relation is 

[0(x),(/.t(x')] ^Vc{^,^'). (77) 

The projector Vc{^, x') will occur frequently in the remainder of this paper, as an expression 
of the fact that the field operator 0(x) and any approximate representations of it must always 
be expressible in terms of the wavefunctions which span the condensate band. 

This projector will play a significant role in the practical implementation of the master 
equations, and is by no means merely a formal requirement. It has two principal effects: 

i) The nonlocality generates a spatial smoothing function, and because the cutoff is in terms 
of energy rather than momentum, the smoothing is stronger where the potential V^(x) is 
larger 

ii) At positions where V{x.) > En, all wavefunctions in the projector are exponentially 
small, so the projector is essentially zero. This means that the boundary conditions on 
any simulational grid are simply that any representation of 0(x) is zero there. 

Although the projector can be written down quite easily using M5\ . this expression is not 
computationally simple — efficient numerical implementation is essential for the practicality 
of any simulations. 

4. Approximate forms of the master equation 

It is not possible to contemplate a numerical solution of the master equation in the form given 
in ( 17 PL but there are two ways in which simplifications can be made, which we will call the 
"quantum optical" master equation and the "high temperature" master equation. 

4.1. The "quantum optical" master equation for a Bose-Einstein condensate 

The method chosen in II15I was to expand the field operators in eigenoperators of Lc- Thus, 
one wrote 

0(x) = ;^X^(x), (78) 

m 

where the eigenoperator requirement is (see fTsl . eq. (28)) 

LcX,„(x) = e™X™(x). (79) 
This method has two disadvantages 
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i) We cannot calculate the operators X,„{x) exactly, although when the Bogoliubov 
approximation is vaUd, they can be expressed in terms of Bogoliubov amplitudes, as 
was done in 1151 . 

ii) The resulting master equation (eq(50a-f) of |l5j) can only be derived by making 
a rotating wave or random phase approximation, whose validity can be questioned. 
However, this master equation is of the Lindblad form (P22J Sect. 5.2.2), and this is 
a highly desirable, though not absolutely essential, property, since it guarantees that 
solutions are density operators with non-negative eigenvalues. 

Useful results have been derived using this form of the master equation. 



4.2. The "high temperature" master equation for a Bose-Einstein condensate 

In this case the essential approximation relies on the fact that eigenfrequencies of Lc are 
small compared to the temperature; precisely 

ri\e.n\lkBT -^l. (80) 

This is a condition which is very often met, to an accuracy of at least 10%, at temperatures 
which one would normally find in a degenerate Rose gases with a significant non-condensate 
fraction. In some sense, the nomenclature "high temperature" is misleading, but in a strict 
sense, it is a fair description of the kind of limit contemplated in (18 0> . The master equation that 
can be derived is not of the Lindblad form, as indeed is also the case for the unapproximated 
master equation ( I70> . However, as shown in 1231 . this probably only means that there are 
transient situations arising after unrealistic (although physically acceptable) initial conditions 
which give rise to a non-positive definite density operator. Such transients usually evolve on 
a time scale more rapid than that used to derive the master equation, so they have no physical 
significance. 

We now develop the approximations based on the condition ( I80> . The functions 
G^^' (u, V, w) and M(u, v,a;) can then be evaluated approximately for fioj <^ ksT, the 
limit in which the stochastic Gross-Pitaevskii equation was originally derived (D- Under this 
condition we can use give 

M(u, V, u;)= (^1 - M(u, V, 0). (81) 

For G(±)(u,v,w) we can take G'+^(u,v,lj) w G^+^(u, v, 0) (Using for example the 
formula (153) of 1151 1 and then write using ( I58> 

G(-) (u, V, c.) « (l - ^) GW (u, V, 0). (82) 

These two expressions will be used to derive all that follows in this paper, but more accurate 
expressions could be used, involving higher powers of huj/kBT, which would extend the 
range of the approximation, as noted by Stoof fT'24' "251. The highest value of e,„ available 
is of the order of Ej^ — /i « 2/i, so that this method would be available for smaller chemical 
potentials and lower temperatures. 

4.2.7. Growth terms Both of G^^^ are sharply peaked as functions of v, so we will also 
approximate (/i(u±v/2) « (j){u.) in the growth terms, and get the approximate master equation 
terms 

Pclgrowth = (l-^) / rf'xG(x){[0(x)pc,0^(x)] ~ Mt(x),0(x)]} 
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+ rf'xG(x) { [(Lc</<(x)) PC, ^^x)] - [pc (ic0t(x)) , 0(x)] } 

+ y d^x G(x) { [0t (x)pc >(x)] ~ [pc</'(x) >^ (x)] } . (83) 
Here we have used the notation 

G(x)= y"d3vG(+)(x,v,0). (84) 

4.2.2. Scattering terms We cannot make the approximation ^(u ± v/2) w 0(u) in 
the scattering terms, because / (Pv M (x, v, 0) is ill defined, since it involves the product 
5{u)i — a;2)(5(Ki — K2). In this case we must keep the full dependence on v: 

PcLcatt = - / d3xd3vAf(x,v,0)|[C/(x + v/2),[C/(x-v/2),pc]] 

' ^ ^C/(x + v/2),[(ict/(x-v/2)),pc]+] 



2kBT 



(85) 



4.2.3. Estimate of the amplitude To understand the nature of the scattering term, one can 
evaluate the Fourier transform 

M(x,k,0)EE^ j d3vc-*-A/(x,v,0), 

^T^TT^ / / d3K2F(Ki,x)[l + i^(K2,x)]J(Ki -K2-k) 

X S{uji — ^2)- (86) 

Since F(K, x) depends on 

eK(x) = ^— + V{^) = hco{K, x), (87) 
2m 

this can be written as 

2u^ f 2m 
M(x,k,0) = ^^-^ y d3KF(eK,x)[l + ^^(eK,x)] — ,5(2K.k-k2). (88) 

We can choose the x-axis in the K-integration parallel to k and and use the delta function to 
eliminate — *■ |k|/2. Using the notation K = {Ky, K^), so that 

we can write 

M(x, k, 0) = j^^^^^^ J d'K fr,^ {e^ik, x)) {l + fr,^ (^^(k, x)) } , (90) 
where 

f^-f^^y^ ^ eiy-,)/LT „ 1. (91) 

Using (|89j, we can change the integration variable to ej^(k, x), where the lower limit of the 
integration is now 

£;^i„(k, x) = max {Er, 1/(x) + n^lkp/Sm) /keT. (92) 
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The integral then can be written 

^''^ {2nfK'\k\ y(iJ.„.„-M)/fc,T [exp(y)-l]2' 

2mu^ 27rmfchT 1 ^^^^ 



(2^)5;i3|j^| ^2 exp((^,„i„(k,x)-/i)/fcBr)-l' 
8a^(fcBT)^ 

(2^)2;i(i?^i„(k,x)-M)|kr ^'^^^ 

The principal dependence on k comes from the l/|k| term, which is most significant for 
small |k|. The operators J7(x ± v/2) are restricted to the condensate band, and thus will 
be significant only for x such that Efi > V{x.). Thus the major contribution comes from 
situations such that £'„iin(k, x) = Ej^, in which case the only k dependence comes from the 
l/|k| term, and we can write 



where we define 



4.2.4. Approximate local form The form \95l is nonlocal, and this may be an important 
feature. Nevertheless, we can try to write an approximate local equivalent in the form (which 
is essentially of the form of the quantum Brownian motion master equation as described in 
l22l Sect. 3.6.1) 



PcLcatt = - / d'xA7(x)|[C/(x),[C/(x),pc]] + ^[C/(x),[(Lcf/(x)),pc] 

(98) 

However, the correct choice of M(x) is not straightforward; a possible estimate can be 
made by taking k = ki — k2, where Tik^ are the momenta of two particles in the 
condensate band. We then compute the average of ( I95> over the region Rc such that 

|ki 2I < \j2m{Eii — y(x)) /h, on the assumption of a noninteracting excitation spectrum: 

' ' ' d^kid^ka, (99) 

(100) 



|ki - k2 



Ra |kl -k2|/ 

6ne{En-vi^)) 



5j2m{ER - y(x)) 



The singularity when Er = V"(x) is harmless, since the occupation is also zero there. Thus 
we can choose 

167ra2(fcBT)2 1 



M(x) 



h{ER-fi) |ki-k2| 



(101) 



4.2.5. Comparison of local and full forms There are obvious technical advantages in using 
the local form ( I98l l of the scattering term instead of the full form (I85> . The local form will 
give the correct stationary distribution for any choice of AI (x), and is in at least this sense 
acceptable. The main difference is that the local form gives essentially the same scattering 
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rate at all momentum transfers, whereas the full form gives a significant drop off at high 
momentum transfers, and this is known to be a significant feature of the kinetics when this is 
treated by forms of the quantum Boltzmann equation, as noted for example in the papers of 
Svistunov, Kagan and Shlyapnikov 1261 1271 . the idea of "flux" in energy space is used, based 
on the locality in energy space of the collision integral. 



5. Wigner function representation 
5.7. Fokker-Planck equations 

Fokker-Planck equations can be derived by using the transforms of |22!!, as described in fP]. 
None of the three principal terms (Hamiltonian, growth, scattering) can be put exactly in the 
form of a genuine probabilistic Fokker-Planck equation, but on the assumption that higher 
order derivatives become less significant, an approximate probabilistic Fokker-Planck, valid 
for large occupations, can be derived. 

The methodology used can be demonstrated for one part of Pclgiowth- 



(x), Pc] - ^ ^Lc(f>{x.)}pc , 0'^(x) 



■ (102) 



The transformation to a Fokker-Planck equation for the Wigner function W of the phase space 
field function a(x) is achieved by the mappings; 

5W 

[Pc,0t(x)] - (103) 
SW 

[0(x),pc] ^ (104) 

^l(t){x)pc Aia(x)l^, (105) 

hLcc^{^)pc ^ Lca[x)W. (106) 

In this equation there are some important points of notation and approximation 

a) The fact that the condensate band is spanned by the finite set of wavefunctions Yn (x) — 
see (I15> — means that we define the phase space amplitude by 

a(x) = ^y„(x)a„, (107) 

n 

(where a„ are independent mode ampUtudes) and that we define a modified functional 
differentiation operator 

This would be a genuine functional differentiation operator if the full set of modes Yn (x) 
were used instead of only the set within the condensate band. Notice that from these 
definitions 

whereas the right hand side is a delta function for a true functional differentiation 
operator. 
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b) Including the effect of the mean field arising from the the non-condensate as an effective 
potential 

Kff(x) = F(x) + 2unNc{^), (110) 
the projected Gross-Pitaevskii operator is introduced as 



Lca{x.) ^Vc \- 



2m 



+ Kff (x)a(x) + u|Q;(x)pa(x] 



(111) 



The projector arises as a consequence of (I103I104I109> etc. 

c) The two replacements (I105I106> are approximate, not only in the sense that higher order 
derivatives are neglected, but also in the sense that all derivatives which should arise 
from the left hand side are neglected. This is not an essential aspect of the method, since 
we can technically still handle the first order derivative terms which would turn up in the 
right hand sides of ( I105ll06f . and these would produce only second order derivatives in 
the final Fokker-Planck equation. This would produce a fractional correction of order 
fi/ksT to the noise terms only. 

The projected Gross-Pitaevskii operator can be written in terms of the Gross-Pitaevskii 
Hamiltonian 



GP 



in the form 



2m 



+ ycff(x)|a(x) 



fl«(x)|^ 



Sa*{x) 

We also define the Gross-Pitaevskii particle number 

|2 



N, 



GP 



J d^x \a{x) 



The term ( I102t then transforms to the term in the Fokker-Planck equation 



d^xG{x)- 



SW 



W Si^lNGP-HGp) 



(112) 



(113) 



(114) 



(115) 



da(x) [Sa*{x) ksT 5a* (x) 

Using the same kind of approximations, the scattering term \9%\ can be transformed to a 
corresponding form, and we obtain, putting all terms together, the Fokker-Planck equation 



dW 



(5a (x) 



d^x- 



+ 



■ 6a* (x) 
d^xd^x M 



iW SHgp 
~YSa*{x) ^ 
iW SHgp 
h (5a (x) 
X + x' 



G(x) 



-G(x) 



SW 
(5a* (x) 
SW 



W SitiNGP-HGp) 



ksT (5a* (x) 
W S{pNgp-Hgp) 



x-x',0 



(5a (x) 
a*(x) 



S 



a*(x') 



SW 



+ 



W SH, 



GP 



(5a* (x') kBTSa*{x') 



— a(x') 



Sa*{x) 
SW 



a(x) 



(5a (x) 
6 



Sa{x) 

W SHgp 

Sa{x') ^ k^5a{x') 



(116) 



The stochastic Gross-Pitaevskii equation: II 



18 



5.2. Stationary solutions of the Fokker-Planck equation 

Independently of the forms of G(x) and A/(x, x'), the stationary solution of the Fokker- 
Planck (II 16> equation is given by 

W^s oc exp ( ^-^^ 1. (117) 

This is the grand canonical distribution expected for a classical field theory whose field a(x, t) 
obeys the Gross-Pitaevskii equation. 

5.3. Stochastic Gross-Pitaevskii equations 



The Fokker-Planck equation (lTT6|l is equivalent to a set of stochastic differential equations, 
which we shall write in the Stratonovich form, in a choice of forms in which various degrees 
of simplification are made of the scattering terms, as follows. The notation (S) to the right of 
an equation indicates that it is in Stratonovich form as in 1281 . 



5.3.7. Full form of the stochastic differential equations Using the full form of the 
scattering term, we have 

(S) (ia(x, t) = — —LcOi{x,t)dt 

+ Vc \ t) - Lcai^, t)} dt + dWa (x, t) 



. J d^x'M (^^^^, X - x', 0^ a(x) 



1_ /■ / x + x' 

2kBT 

X |a*(x',i) (Zca(x',t)) -a(x',t) (Zca(x', dt 



+ ia(x)dW^M(x,t)|. (118) 

The noise dWoi^, ^) is complex, while dWmi^, t) is real; they are independent of each other, 
and satisfy the relations 

dW^{x,t)dWG{x',t) = 2G'(x)(5(x - x') dt, (119) 
dVFG(x,t)dM^G(x',t) = dWa{x,t)dWG{x,t) =0, (120) 

dWM (x, t) dWM (x', t) = 2M ( ^^4r-. X - x', ) dt. (121) 



5.3.2. Simplified non-local form The implementation of the last two lines of ( II 18l l 
obviously presents some technical issues. However, if we simplify ( 195 > by setting 
Si„in(k, x) Efi (as discussed there), the nonlocal form can be considerably simplified. 
The operator 1 /\/— V^, of which 1/ 1 k| is the Fourier transform, is not singular when acting on 
functions with well behaved Fourier transform as k ^ 0. Using this, the stochastic differential 
equation becomes 

(S) da(x,t) = — —LcOi(x,t)dt 
n 

+ Vcl {Ma(x,t) - hLca{x,t)} dt + dWG{^,t) 
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- ^1^^^^ {a*(x,i)Zca(x,t) - a(x,0 (Zca(x,t))*} dt 

+ ia(x)dW^^(x,t)|. (122) 

The noise dWc (x, t) is complex, while dWM (x, t) is real; they are independent of each other, 
and satisfy the relations 

dW^{x,t)dWG{yi' ,t) = 2G'(x)(5(x-x')dt, (123) 
dWG{-^,t)dWG{^ ,t) = dVK5(x,<)dVFG(x',t) = 0, (124) 

(iWA4(x,t)dt¥M(x',t) 6{x-x')dt. (125) 

v— 

5.5.5. Local form of the stochastic differential equations These take the form 

(S) da{-K,t) ^ - ^Lc(a(x,i)) dt 

+ 'Pc\ - hLca{x,t)} dt + dWci^^t) 

M(x)a(x) r», ,- , , , , f~ , 
IkeT \ i)^ca(x, t) - a(x, t) (Lca(x, i)) | 

+ ia(x)dW^^7(x,t)|. (126) 

The noise dWciyi-^ t) is complex, while (iW^7(x, t) is real; they are independent of each other, 
and satisfy the relations 

dVF5(x,i)dW^G(x',i) = 2G(x),5(x-x')di, (127) 
dW^G(x,t)dW^G(x',i) dVKG(x,t)dW^5(x',t) = 0, (128) 
dM^£f(x, dM^A7(x', t) = 2A7(x) ,5(x - x') dt. (129) 



5.4. Comparison with other methods 

The stochastic differential equations we have derived have similarities with those we have 
previously derived, as well as with Stoof's. However, our way of implementing the idea of 
eliminating higher energy thermalized modes has significant differences from Stoof's: 

i) There are major technical differences in how the elimination is done. Ours is based on the 
ideas of quantum optics, which are used to develop a quantum mechanical master equation. 
The master equation is then transformed using the Wigner function to give a Fokker-Planck 
equation, which is equivalent to a stochastic differential equation. Stoof uses a functional 
integral formulation of the Keldysh method, in which the elimination is achieved in the action 
integral. This method is almost certainly equivalent to our quantum optical method. Thus, 
although the technical methods used appear very different, this is not physically significant. 

ii) However the choice of what to eliminate is different. Stoof eliminates modes with energy 
> /i, and finds a Fokker-Planck equation (equivalent to a stochastic differential equation) 
which involves self energy functions, which he evaluates using the many-body T-matrix 
method. In contrast we carry out the elimination as a two-stage process, as detailed in 
Sects. 12131 Thus we separate the elimination process which is required to give an effective 
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field theory from that required to give a master equation for the condensate band. This has the 
advantage that we do not need to use a many body T-matrix description. 

iii) The equations which arise are given above as ( I118I122I126> . depending on degrees of 
approximation used. The principal differences are the appearance of the projector, and the 
inclusion of scattering terms. There is also the implicit difference that the field variable we 
use is not exactly the same as Stoof 's, since it includes a wider range of states. 

iv) Our earlier work [IJ included the scattering terms, but was unable to give any precise 
value to them because of the crudeness of the methods used, and also for reasons related to 
the fact, which we have found in this paper, that the scattering cannot be described locally 
without making very crude approximations, such as those which lead to the form ( I126> . The 
explicit appearance of the projector in the equations did not occur there or in Stoof 's method. 
It is possible that in numerical implementation, the projector will not always be essential, but 
we are confident that there are situations in which its inclusion will be important. 

6. Conclusions and outlook 

The problem we have set ourselves in this paper is to produce a description of finite 
temperature condensed or nearly condensed Bose gases which is both accurate and 
implementable. Apart from those used to derive the basic Markovian master equation ( I33> . 
there are only two significant approximations made in our derivation: 

i) The energy eigenvalues of the excitations are small compared to temperature, as noted in 
Sect. l4.2l This is needed to get a master equation of a reasonably simple kind, the "high 
temperature master equation." 

ii) The occupations of the modes being treated in the condensate band are significantly 
larger than unity, as noted in Sect. l5.ll This is essential to be able to use a Wigner 
function representation of the master equation which becomes equivalent to a classical 
field representation. 

These two conditions are essentially the same, since the occupation of a mode of energy 
e is large if e <C fc^T, even though the reasons for the conditions are quite independent. 
Therefore, it is clear that if one is using the "high temperature master equation," then it is also 
sensible to use the Wigner function classical field representation, since the two are accurate 
to the same degree of approximation. 

6.1. The validity of the classical field representation 

Classical field representations have often been introduced heuristically by the argument that 
the quantum field operator can be replaced by a classic field function provided the occupations 
of the modes are high. The Wigner function methodology which we use gives systematic 
way of implementing these heuristic ideas, and gives a way a assessing the validity of the 
formalism^. It is important to realize that, when formulated through the Wigner function, the 
classical field method gives a truly quantum mechanical description of the system, subject 
only to the technical approximations made. This means that we expect the major quantum 
mechanical aspects to be correctly treated; in other words, the classical field method, correctly 
and carefully formulated, has no heuristic approximations introduced to provide refuge from 

I The alternative P- and Q-function descriptions are better called phase space representations, since the resulting 
equations of motion possess non-negligible noise terms of a purely quantum nature, and thus the field cannot be said 
to be a classical field. 
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quantum mechanics — it is a valid way of treating quantum mechanics in a certain degree of 
approximation, which fortunately behaves very much like a classical description. 

As a result, certain quantum properties make their presence felt when using this method. 
The Wigner function exhibits "vacuum occupation" of the modes — the mean square of the 
classical field is non-zero even when there are no particles. This follows from the relationship 
between the classical field averages and the quantum averages in the form 



which says that the minimum occupation exhibited by the classical field is the "vacuum 
occupation," half a particle per mode — an expression of the Heisenberg uncertainty principle. 
When summed over the infinite number of modes which constitute the full system, this gives 
a divergent field function; thus a cutoff must be introduced, such as the one we have chosen, 
which defines the boundary of the condensate band. 

There are two different ways in which the classical field method applied to the condensate 
band can be valid. 

i) The temperature may be so low that there is little occupation at the top of the condensate 
band, so that all of the noise terms derived in Sect.|5]are zero. At low energies, the 
excitations of a Bose condensate are largely phonon-like, and amount to quantized 
shape and density oscillations of the condensate itself. These are only noticeable if 
their occupation is large; modes with low occupation, although badly described by the 
classical field method, provide a completely negligible contribution to the system. 
However, since there is no external noise, in this case the initial occupations of the 
modes must be chosen so as to represent the existence of a finite temperature. Even 
if the temperature is zero, these occupations must be chosen to give the correct "vacuum 
occupation," and this gives rise to irreversible effects, as first noted by Steel et al. fZ9] . 

ii) If the temperature is not very low, the total number of particle in modes with occupations 
< 1 can be very significant, and this is necessarily the case during the process of 
condensate growth from a vapour. The influence of these modes has to be included, 
and this is what has been done in this paper. 

In principle the initial conditions should be chosen as in i) to give the right "vacuum 
occupation", but their effect rapidly dies out because of the irreversible coupling to the 
reservoir which generates the noise and damping terms. 

6.2. The degenerate non-condensed Bose gas 

It is known that during the process of condensate growth the occupations of a large number 
of modes become significantly larger than one before the appearance of a single dominant 
mode, the condensate. These modes are of course describable by our stochastic classical 
field equations (I118H129> , since we have nowhere assumed the existence of a condensate, and 
because their high occupation makes the classical field method valid. Kagan, Svistunov and 
Shlyapnikov 1271 in their pioneering work on the initiation of Bose-Einstein condensation, 
also argued that a description in terms of the Gross-Pitaevskii equation was possible for such 
a system. However, because their philosophy was more qualitative, the technical details which 
we are compelled to attend to do not appear in their work. 

6.3. The projector into the condensate band 

One of the essential features of the description is the presence of a projector in the equations 
of motion for the classical field. This prevents non-condensate band modes being wrongly 




(130) 
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included in numerical simulations. Such a projector has already been implemented in 
dJ |9] 1101 . but only for a homogeneous system. While the projector for a system with a 
trap potential can be written down explicitly as in ( I15I16> . it is not an easy task to implement 
this projector efficiently. 

The projector deals correctly with two problems which arise in practical simulations: 

i) The fineness of the mesh: The distance between mesh points determines the highest 
momentum which can be represented in the simulation. However, one cannot choose this 
to be the physical cutoff, even in the case of no trapping potential, since nonlinear terms 
such as Vc {|q^(x)Pq:(x)} in il 1 1> are wrongly represented this way. To represent this 
term correctly without aliasing, the numerical wavenumber cutoff must extend to at least 
two times the physical cutoff. If one assumes the physical cutoff is to be given by the 
mesh fineness, incorrect results will be obtained unless there is no significant occupation 
above one half of the mesh wavenumber cutoff. 

In addition, when a trap potential is present, a cut at a minimum wavelength does not 
correspond to a cut at definite energy. This would mean that even non-interacting atoms 
would pass from the condensate band to the non-condensate band with the progression 
of time, presenting some difficulties for the formalism 

ii) The boundary condition at the edge of the spatial grid: Because the projector involves 
trap eigenfunctions only up to energy Ej^, in the case of a non-vanishing trap potential 
there is a distance R such that for |x| > R projected functions become exponentially 
small. This means that not only are the field functions ct{x.) of Sect.|5] exponentially 
small there, but also the added noise terms as well. Therefore in practical simulations 
on a grid of finite size, the issue of the boundary condition for the noise and the field 
function at the edge of the grid is determined; the simulation region must be so large that 
all field and noise functions can be set equal to zero at the boundary. 

Where there is no trap potential, the box inside which the simulation is being 
implemented has real physical significance, and periodic boundary conditions are usually 
chosen. 

6.4. Applications 

i) Condensate growth: The theory of condensate growth is at this stage not entirely 
satisfactory — there is still disagreement between theory and experiment under certain 
conditions 1301 1311 . The main defect in computations of condensate growth 1301 1311 
l32l l33l l34l l35l 1361 to date is an inadequate treatment of the phonon-like quasiparticle 
excitations. Using the formulation proposed here would definitely include these 
correctly. 

ii) Vortex lattice growth: The theoretical situation in the growth of vortex lattices is in 
a very preliminary stage. The stabilization of the vortices into a lattice is clearly a 
result of some kind of irreversible process, for which a phenomenological model was 
first proposed by Tsubota et al. 1371 . while we presented a physically based model 1381 
based on a rotating frame version of our phenomenological growth equation (l] which 
attributed the necessary irreversibility to interaction with a thermal cloud. 

In contrast, Lobo et al. i39j used a simple stochastic classical field model of the type 
mentioned above, in which the noise arises from initial conditions with no interaction 
with a thermal cloud to provide the requisite damping. However, the noise in their 
implementation of this model is nonzero throughout the simulation region, and this 
causes some difficulty in deciding the appropriate boundary condition at the edge of the 
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simulation region, since they choose the high energy cutoff to be given by the fineness 
of the simulation mesh, with no projector as in this paper In effect, this means that 
the vacuum can transfer angular momentum to the system, and that indeed the vacuum 
is different depending on whether the simulation is carried out in a rotating frame or a 
non-rotating frame. 

Our proposed methodology would avoid the boundary condition difficulties of the last 
model, combining its ideas with those of our earlier model. 

ii) Heating of a condensate by mechanical disturbance: When a condensate is 
mechanically disturbed, some heating occurs. This is clearly also a phenomenon which 
can be treated by our methodology. 

6.5. Technical aspects: 

The proposed stochastic differential equations J118H1291 are not as simple as one would 
like, but are nevertheless probably quite practical even in three dimensions. The noise and 
damping arise due to the growth and scattering terms in the master equation. The former 
describe the transfer of particles between the condensate and non-condensate bands, while 
the latter represents the effect of non-condensate particles colliding with condensate band 
atoms with one particle remaining in each band. The growth terms in the stochastic Gross- 
Pitaevskii equation are relatively easy to implement numerically, but there are definitely some 
challenges for the scattering terms. The essential feature of the projector into the condensate 
band in the equations is an additional challenge, relatively easy to conceive, but formidable to 
implement. However, the projector is not of merely cosmetic significance, and this challenge 
must be faced. 

6.6. Outlook 

The practical implementation of our methodology will be the feature of forthcoming papers 
which will include treatments of condensate growth, vortex lattices and heating. 
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